open(FH,"/home/uec-00/yapingli/software/gatk_read/Sting/resources/Normal_SNP.filterCNV.txt") or die;
my @in=<FH>;
chomp(@in);
close(FH);
my $count = 0;
my $count_het = 0;
my $count_array_het = 0;
my $count_hom = 0;
my $fp = 0;
my $fn = 0;
my $total_count=0;
my $ag_count=0;
my $count_array_hom = 0;
open(FH,"<normalMerge.m-15.b-1.bs0.04.bisulfiteRate-0.97.performace.onlyvariant.raw.vcf") or die;
my @vcf=<FH>;
chomp(@vcf);
close(FH);
my %vcf_hash = ();
foreach my $line(@vcf){
	next if $line=~/^\#/;
	my @splitvcf=split "\t",$line;
#	next if $splitvcf[5] <=30;
	$vcf_hash{$splitvcf[2]} = $line;
}
foreach my $line(@in){
	my @splitin=split "\t",$line;
	next if(($splitin[5] eq "A" && $splitin[6] eq "A") || ($splitin[5] eq "A" && $splitin[6] eq "T") || ($splitin[5] eq "T" && $splitin[6] eq "A") || ($splitin[5] eq "T" && $splitin[6] eq "T"));
	if($splitin[3] ne $splitin[4]){
		$count_array_het++;
	}
	else{
		$count_array_hom++;
	}
	if( exists $vcf_hash{$splitin[0]} ){
		$count++;
		my @splitvcf=split "\t",$vcf_hash{$splitin[0]};
		if($splitvcf[9]=~/0\/1/){
			$count_het++;
			if($splitin[3] eq $splitin[4]){
				$fp++;
			#	print("fp:\t$splitin[0]\t$splitin[1]\t$splitin[2]\n");
			}
		}
		else{
			if($splitvcf[4]=~/[AGCT]\,/){
				$count_het++;
				if($splitin[3] eq $splitin[4]){
					$fp++;
				#	print("fp:\t$splitin[0]\t$splitin[1]\t$splitin[2]\n");
				}
			}
			else{
				if($splitin[3] ne $splitin[4]){
					$fn++;
					if($line=~/\d+\tA\tB\tA\tG/ or $line=~/\d+\tA\tB\tG\tA/ or $line=~/\d+\tA\tB\tT\tC/ or $line=~/\d+\tA\tB\tC\tT/){
						$ag_count++;
					}
				#	print("false neg:\t$splitin[0]\t$splitin[1]\t$splitin[2]\n");
				}
				$count_hom++;
			}
			
		}
	}
	else{
		if($splitin[3] ne $splitin[4]){
			$fn++;
			if($line=~/\d+\tA\tB\tA\tG/ or $line=~/\d+\tA\tB\tG\tA/ or $line=~/\d+\tA\tB\tT\tC/ or $line=~/\d+\tA\tB\tC\tT/){
				$ag_count++;
			}
			#print("false neg:\t$splitin[0]\t$splitin[1]\t$splitin[2]\n");
		}
		$count_hom++;
		##print("$splitin[2]\n");
	}
	$total_count++;
}
my $fp_rate = $fp/$count_array_hom;
my $fn_rate = $fn/$count_array_het;
print("count is: $count\n");
print("total count is: $total_count\n");
print("het count is: $count_het\n");
print("array het count is: $count_array_het\n");
print("array hom count is: $count_array_hom\n");
print("false pos num is: $fp\t fp rate is: $fp_rate\n");
print("false neg num is: $fn\t fn rate is: $fn_rate\n");
print("ag count in fn is: $ag_count\n");

open(FH,"/home/uec-00/yapingli/code/performance_test/k-allele/performance-normal-15-1-1-1-0.1.o10818201") or die;
my @kAllele=<FH>;
chomp(@kAllele);
close(FH);

my $fn_kallele=0;
my $fp_kallele=0;
foreach my $line(@kAllele){
	next if($line !~ /\S+\t\S+\t\S+\t\S+\t\S+\t\S+\t\S+/);
	my @splitin=split "\t",$line;
	next if(($splitin[5] eq "A" && $splitin[6] eq "A") || ($splitin[5] eq "A" && $splitin[6] eq "T") || ($splitin[5] eq "T" && $splitin[6] eq "A") || ($splitin[5] eq "T" && $splitin[6] eq "T"));
	if(($splitin[5] eq "G" && $splitin[6] eq "G") || ($splitin[5] eq "C" && $splitin[6] eq "C")){
		$fp_kallele++;
	}
	else{
		$fn_kallele++;
	}
	
}
my $fp_rate_kallele = $fp_kallele/$count_array_hom;
my $fn_rate_kallele = $fn_kallele/$count_array_het;
print("false pos num in k allele method is: $fp_kallele\t fp rate is: $fp_rate_kallele\n");
print("false neg num in k allele method is: $fn_kallele\t fn rate is: $fn_rate_kallele\n");
